Rachel_QTL_GWAS
Rachel Olson
March 7, 2020
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.2.1 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.4
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: lattice## Loading required package: crayon##
## Attaching package: 'crayon'## The following object is masked from 'package:ggplot2':
##
## %+%## Loading required package: Biobase## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following object is masked from 'package:Matrix':
##
## which## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Parsed with column specification:
## cols(
## Chr = col_character(),
## Position = col_double(),
## Marker = col_character(),
## ref = col_character(),
## alt = col_character()
## )Manipulate Files
Read the Cross
## --Read the following data:
## 571 individuals
## 30328 markers
## 10 phenotypes## Warning in summary.cross(cross): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
## (Perhaps it is in basepairs?)## --Cross type: f2## Warning in summary.cross(crossobj): Some chromosomes > 1000 cM in length; there may be a problem with the genetic map.
## (Perhaps it is in basepairs?)## F2 intercross
##
## No. individuals: 571
##
## No. phenotypes: 10
## Percent phenotyped: 100 94.7 94.7 87 99.8 99.8 99.8 96.5 96.5 96.5
##
## No. chromosomes: 7
## Autosomes: Fvb1 Fvb2 Fvb3 Fvb4 Fvb5 Fvb6 Fvb7
##
## Total markers: 30328
## No. markers: 3890 4101 4745 4013 4451 5381 3747
## Percent genotyped: 100
## Genotypes (%): AA:29.2 AB:35.6 BB:35.2 not BB:0.0 not AA:0.0Optional: form linkage groups #{r, echo=F} newcrossobj <- formLinkageGroups(crossobj, max.rf=0.35, min.lod=6, reorgMarkers=TRUE, verbose=FALSE) newcrossobj <- orderMarkers(newcrossobj, use.ripple = TRUE, window = 5, map.function = c("kosambi"), verbose = T) plot(newcrossobj) summaryMap(newcrossobj) #
QTL Mapping and Linkage
## [1] 0.4960421## [1] 0.3184846 0.9235690
## ind1 ind2 prop_match
## 1 52QKSimTco R9KUKOjnBl 0.924## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 30 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## chr pos CWT_WEO CWT_SAL1 CWT_SAL2
## AX-184840405 Fvb1 12035 26.0 53.8 60.2
## AX-184687419 Fvb2 18860 25.8 46.1 52.7
## AX-184908616 Fvb3 22445 26.7 47.9 56.4
## AX-184169910 Fvb4 3765 27.8 46.7 54.5
## AX-184843747 Fvb5 295 27.4 55.2 65.8
## AX-184677344 Fvb6 2525 27.5 52.3 60.3
## AX-184202548 Fvb7 4820 28.1 57.4 65.6## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 30 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## chr pos CCT_WEO CCT_SAL1 CCT_SAL2
## AX-184346603 Fvb1 8420 13.0 14.4 16.48
## AX-184687419 Fvb2 18860 14.2 16.5 18.85
## AX-184908616 Fvb3 22445 11.0 10.1 13.93
## AX-184169910 Fvb4 3765 11.5 12.7 15.78
## AX-184170905 Fvb5 6970 11.4 12.6 11.26
## AX-89849903 Fvb6 10950 11.1 9.7 7.44
## AX-184164903 Fvb7 10865 11.9 13.1 12.10## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 74 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[1], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 1 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.## Warning in checkcovar(cross, pheno.col, addcovar, intcovar, perm.strata, : Dropping 20 individuals with missing phenotypes.## Warning in scanone(cross, chr, pheno.col[i], model, method, addcovar,
## intcovar, : First running calc.genoprob.
## chr pos CAWT_WEO CAWT_SAL1 CAWT_SAL2
## AX-184335455 Fvb1 14285 31.3 70.0 70.9
## AX-166511927 Fvb2 3805 28.6 84.9 92.9
## AX-184612678 Fvb3 9615 29.6 90.8 98.2
## AX-123367232 Fvb4 19085 29.4 50.6 48.3
## AX-184595061 Fvb5 20095 29.4 66.8 67.7
## AX-184615096 Fvb6 23015 30.2 86.9 92.3
## AX-184590330 Fvb7 180 27.1 81.7 96.9GWAS
## [1] 1.02
WEO
## iteration LogLik wall cpu(sec) restrained
## 1 -188.263 6:47:24 1 0
## 2 -175.97 6:47:24 1 0
## 3 -173.329 6:47:25 2 0
## 4 -172.915 6:47:26 3 0
## 5 -172.892 6:47:27 4 0
## 6 -172.89 6:47:28 5 0
## 7 -172.89 6:47:28 5 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -222.416 6:50:28 1 0
## 2 -214.567 6:50:28 1 0
## 3 -211.774 6:50:29 2 0
## 4 -211.138 6:50:30 3 0
## 5 -211.076 6:50:30 3 0
## 6 -211.069 6:50:31 4 0
## 7 -211.069 6:50:31 4 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -164.067 6:53:21 1 0
## 2 -158.203 6:53:21 1 0
## 3 -157.561 6:53:22 2 0
## 4 -157.52 6:53:22 2 0
## 5 -157.519 6:53:23 3 0
## Performing GWAS evaluation
SAL1
## iteration LogLik wall cpu(sec) restrained
## 1 -138.928 6:55:49 0 0
## 2 -115.797 6:55:50 1 0
## 3 -114.631 6:55:51 2 0
## 4 -114.573 6:55:52 3 0
## 5 -114.573 6:55:53 4 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -213.835 6:58:32 0 0
## 2 -198.408 6:58:33 1 0
## 3 -193.384 6:58:34 2 0
## 4 -191.709 6:58:35 3 0
## 5 -191.323 6:58:35 3 0
## 6 -191.225 6:58:36 4 0
## 7 -191.2 6:58:37 5 0
## 8 -191.194 6:58:37 5 0
## 9 -191.192 6:58:38 6 0
## 10 -191.192 6:58:39 7 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -54.232 7:1:19 1 0
## 2 -26.0886 7:1:20 2 0
## 3 -24.3347 7:1:21 3 0
## 4 -24.087 7:1:21 3 0
## 5 -24.073 7:1:22 4 0
## 6 -24.0721 7:1:23 5 0
## Performing GWAS evaluation
SAL2
## iteration LogLik wall cpu(sec) restrained
## 1 -110.883 7:4:3 0 0
## 2 -86.6082 7:4:4 1 0
## 3 -85.8537 7:4:5 2 0
## 4 -85.8348 7:4:6 3 0
## 5 -85.8347 7:4:6 3 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -208.738 7:6:42 1 0
## 2 -193.975 7:6:42 1 0
## 3 -188.208 7:6:43 2 0
## 4 -186.261 7:6:43 2 0
## 5 -185.89 7:6:44 3 0
## 6 -185.822 7:6:45 4 0
## 7 -185.809 7:6:45 4 0
## 8 -185.807 7:6:46 5 0
## 9 -185.807 7:6:47 6 0
## Performing GWAS evaluation
## iteration LogLik wall cpu(sec) restrained
## 1 -25.7674 7:9:41 1 0
## 2 6.386 7:9:42 2 0
## 3 7.01255 7:9:42 2 0
## 4 7.06207 7:9:43 3 0
## 5 7.06264 7:9:44 4 0
## Performing GWAS evaluation
Checking the outcome
## [1] 1.001463
## [1] 1.037738
## [1] 1.095772
## [1] 1.04106
## [1] 0.9822432
## [1] 0.9875107
## [1] 1.177989
## [1] 1.096744
## [1] 1.086727
Mutivariate attempt #```{r, echo=FALSE} #CWT_multiple years miss10 <- pheno[is.na(pheno$CWT_SAL2),] MixCWT_yr <- sommer::GWAS(fixed = cbind(CWT_SAL1,CWT_SAL2)Year, random=vs(ID, Gu= KA), rcov=~units, data=pheno[which(pheno$ID %in% rownames(KA)),], n.PC = 0, min.MAF=0, M= GO[!(rownames(GO) %in% miss10$ID),], gTerm = “u:ID”)
us7 <- as.data.frame(t(UniCWT_SAL2\(scores)) us7\)Marker <- rownames(us7) MI2_CWT_SAL2 <- merge(markers,us7,by=“Marker”,all.x = TRUE) colnames(MI2_CWT_SAL2) <- c(“Marker”,“Chrom”,“Position”,“ref”,“alt”,“CWT_SAL2 beta”,“p.val”,“CWT_SAL2 Fstat”,“R2”,“R2s”) manhattan(MI2_CWT_SAL2, pch=20,cex=.5, PVCN = “color score”) ```